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Abstract 

In a previous study [DD07b] we presented a novel visco-potential free surface flows 
formulation. The governing equations contain local and nonlocal dissipative terms. 
From physical point of view, local dissipation terms come from molecular viscos- 
ity but in practical computations, rather eddy viscosity should be used. On the 
other hand, nonlocal dissipative term represents a correction due to the presence 
of a bottom boundary layer. Using the standard procedure of Boussinesq equations 
derivation, we come to nonlocal long wave equations. In this article we analyse 
dispersion relation properties of proposed models. The effect of nonlocal term on 
solitary and linear progressive waves attenuation is investigated. Finally, we present 
some computations with viscous Boussinesq equations solved by a Fourier type 
spectral method. 
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1 Introduction 



Even though the irrotational theory of free-surface flows can predict success- 
fully many observed wave phenomena, viscous effects cannot be neglected 
under certain circumstances. Indeed the question of dissipation in potential 
flows of fluid with a free surface is an important one. As stated by [LH92], 
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it would be convenient to have equations and boundary conditions of com- 
parable simplicity as for undamped free-surface flows. The peculiarity here 
lies in the fact that the viscous term in the Navier-Stokes (NS) equations is 
identically equal to zero for a velocity deriving from a potential. There is also 
a problem with boundary conditions. It is well known that the usual non-slip 
condition on the solid boundaries does not allow to simulate free surface flows 
in a finite container. Hence, some further modifications are required to permit 
the free surface particles to slide along the solid boundary. These difficulties 
were overcome in our recent work [DD07b] and [DDZ08]. 

The effects of viscosity on gravity waves have been addressed since the end 
of the nineteenth century in the context of the linearized Navier-Stokes (NS) 
equations. It is well-known that Lamb [Lam32] studied this question in the 
case of oscillatory waves on deep water. What is less known is that Boussinesq 
studied this effect as well [Bou95]. In this particular case they both showed 
that 

— = -2uk^a, 1 
dt ^ ^ 

where a denotes the wave amplitude, u the kinematic viscosity of the fluid 

and k = 2tt / 1 the wavenumber of the decaying wave. Here £ stands for the 

wavelength. This equation leads to the classical law for viscous decay, namely 



a{t) = aoe 



-2uk'^t 



Let us consider a simple numerical application with g = 9.8 m/s"^, i = 3 m 
and molecular viscosity u = 10~® m?/s. According to the formula (2), this 
wave will take to ~ 8 x 10^ s or about one day before losing one half of its 
amplitude. This wave will attain velocity 

'— ^ 2.16 ml s 

k 

and travel the distance equal to L = Cp ■ to ~ 170 km. This estimation is 
exaggerated since the classical result of Boussinesq and Lamb does not take 
into account energy dissipation in the bottom boundary layer. We will discuss 
the question of linear progressive waves attenuation in Section 5. Another 
point is that the molecular viscosity u should be replaced by eddy viscosity 
ft which is more appropriate in most practical situations (see Remark 6 for 
more details). 

The importance of viscous effects for water waves has been observed in various 
experimental studies. For example, in [ZG71] one can read 

. . . However, the amplitude disagrees somewhat, and we suppose that this 
might be due to the viscous dissipation. . . 

[WuSl] also mentions this drawback of the classical water wave theory: 
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. . . the peak amplitudes observed in the experiments are shghtly smaller 
than those predicted by the theory. This discrepancy can be ascribed to the 
neglect of the viscous effects in the theory. . . 

Another example is the conclusion of [BPS81]: 

... it was found that the inclusion of a dissipative term was much more im- 
portant than the inclusion of the nonlinear term, although the inclusion of 
the nonlinear term was undoubtedly beneficial in describing the observa- 
tions. . . 

Another source of dissipation is due to bottom friction. An accurate computa- 
tion of the bottom shear stress f is crucial for calculating sediment transport 
fluxes. Consequently, the predicted morphological changes will greatly depend 
on the chosen shear stress model. Traditionally, this quantity is modelled by 
a Chezy-type law 

Az=-h = Cfp\Uh\Uh, 

where Uh = u{x, y, —h, t) is the fluid velocity at the bottom, C/ is the friction 
coefficient and p is the fluid density. Two other often used laws can be found 
in [DD07a, Section 4.3], for example. One problem with this model is that r 
and Uh are in phase. It is well known [LSVO06] that in the case of a laminar 
boundary layer, the bottom stress is | out of phase with respect to the 

bottom velocity. 

Water wave energy can be dissipated by different physical mechanisms. The 
research community agrees at least on one point: the molecular viscosity is 
unimportant. Now let us discuss more debatable statements. For example if we 
take a tsunami wave and estimate its Reynolds number, we find Re ~ 10^. So, 
the flow is clearly turbulent and in practice it can be modelled by various eddy 
viscosity models. On the other hand, in laboratory experiments the Reynolds 
number is much more moderate and sometimes we can neglect this effect. 
When nonbreaking waves feel the bottom, the most efficient mechanism of 
energy dissipation is the bottom boundary layer. This is the focus of our paper. 
We briefly discuss the free surface boundary layer and explain why we do not 
take it into account in this study. Finally, the most important (and the most 
challenging) mechanism of energy dissipation is wave breaking. This process is 
extremely difficult from the mathematical but also the physical and numerical 
points of view since we have to deal with multivalued functions, topological 
changes in the flow and complex turbulent mixing processes. Nowadays the 
practitioners can only be happy to model roughly this process by adding ad- 
hoc dissipative terms when the wave becomes steep enough. 

In this work we keep the features of undamped free-surface flows while adding 
dissipative effects. The classical theory of viscous potential flows is based on 
pressure and boundary conditions corrections [JW04] due to the presence of 
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viscous stresses. We present here another approach. 

Currently, potential flows with ad-hoc dissipative terms are used for example 
in direct numerical simulations of weak turbulence of gravity waves [DKZ03, 
DKZ04, ZKPD05]. There have also been several attempts to introduce dis- 
sipative effects into long wave modelling [Mei94, DD07a, CG07, Kha97]. We 
would like to underline that the last paper [Kha97] contains also a nonlocal 
dissipative term in time. 

The present work is a direct continuation of the recent studies [DDZ08, DD07b] . 
In [DDZ08] the authors considered periodic waves in infinite depth and deriva- 
tion was done in two-dimensional (2D) case, while in [DD07b] we removed 
these two hypotheses and all the computations are done in 3D. This point is 
important since the vorticity structure is more complicated in 3D. In other 
words we considered a general wavetrain on the free surface of a fluid layer of 
finite depth. As a result we obtained a formulation which contains a nonlocal 
term in the bottom kinematic condition. Later we discovered that the nonlocal 
term in exactly the same form was derived in [LO04]. The inclusion of this 
term is natural since it represents the correction to potential flow due to the 
presence of a boundary layer. Moreover, this term is predominant since its 
magnitude scales with 0{^\fv\ while other terms in the free-surface boundary 
conditions are of order 0{y). The importance of this effect was pointed out in 
the classical literature on the subject [Lig78]: 

. . . Bottom friction is the most important wherever the water depth is sub- 
stantially less than a wavelength so that the waves induce significant hor- 
izontal motions near the bottom; the associated energy dissipation takes 
place in a boundary layer between them and the solid bottom. . . 

This quotation means that this type of phenomenon is particularly important 
for shallow water waves like tsunamis, for example [DD06, Dut07]. Here we 
present several numerical computations based on the newly derived governing 
equations and analyse dispersion relation properties. 

We would like to mention here a paper of N. Sugimoto [SugQl]. The author 
considered initial-value problems for the Burgers equation with the inclusion 
of a hereditary integral known as the fractional derivative of order |. The 
form of this term was justified in previous works [Sug89, SugQO]. Note, that 
from fractional calculus point of view our nonlocal term (6) is also a half-order 
integral. 

Other researchers have obtained nonlocal corrections but they differ from ours 
[KM75]. This discrepancy can be explained by a different scaling chosen by 
Kakutani & Matsuuchi in the boundary layer. Consequently, their governing 
equations contain a nonlocal term in space. The performance of the present 
nonlocal term (6) was studied in [LSVO06]. The authors carried out in a wave 
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tank a set of experiments, analyzing the damping and shoaling of solitary 
waves. It is shown that the viscous damping due to the bottom boundary 
layer is well represented by the nonlocal term. Their numerical results fit well 
with the experiments. The model not only properly predicts the wave height 
at a given point but also provides a good representation of the changes on 
the shape and celerity of the soliton. We can conclude that the experimental 
study by P. Liu et al. [LSVO06] validates this theory. 

The present article is organized as follows. In Section 2 we estimate the rate of 
viscous dissipation in different regions of the fluid domain. Then, we present 
basic ideas of derivation and come up with visco-potential free-surface flows 
formulation. At the end of Section 3 we give corresponding long wave models: 
nonlocal Boussinesq and KdV equations. Section 4 is completely devoted to 
the analysis of linear dispersion relation of complete and long wave models 
introduced in previous section. Last two sections deal with linear progressive 
and solitary waves attenuation respectively. Finally, the paper is ended by 
some conclusions and perspectives. 



2 Anatomy of dissipation 

In this section we briefly discuss the contribution of different flow regions into 
water wave energy dissipation. We conventionally [Mei94] divide the flow into 
three regions illustrated on Figure 1. On this figure Sf and 5";, stand for free 
surface and bottom respectively. Then, Rf and Rb denote the interior 
region, free surface and bottom boundary layers. 




Figure 1. Conventional partition of the flow region into interior region and free 
surface, bottom boundary layers. 

In order to make some estimates we introduce the notation which will be used 
in this section: fj, is the dynamic viscosity, 6 = 0[^/Ji) is the boundary layer 
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thickness, to is the characteristic time, oq is the characteristic wave amphtude 
and i is the wavelength. 

We assume that the flow is governed by the incompressible Navier-Stokes 
equations: 



V-n = 0, 

— + u ■ Vn + -Vp = 9 + -V -T, 
at p P ~ 

where r is the viscous stress tensor 

1 f dui duj 



We multiply the second equation by u and integrate over the domain Vt with 
boundary dVt to get the following energy balance equation: 



n an 

= J + T^n ■ u da + J pg ■ u dVt — — J t: t dil . 

an n q 

V ' 

r 

In this identity each term has a precise physical meaning. The left-hand side 
is the total rate of energy change in Q. The second term is the flux of energy 
across the boundary. On the right-hand side, the first integral represents the 
rate of work by surface stresses acting on the boundary. The second integral 
is the rate of work done by the gravity force throughout the volume, and the 
third integral T is the rate of viscous dissipation. We focus our attention on 
the last term T. We estimate the order of magnitude of the rate of dissipation 
in various regions of the fluid. 

We start by the interior region Ri. Outside the boundary layers, it is reasonable 
to expect that the rate of strain is dominated by the irrotational part of the 
velocity whose scale is ^ and the length scale is the wavelength i. The energy 
dissipation rate is then 



Inside the bottom boundary layer the normal gradient of the solenoidal part 
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of u dominates the strain rate, so that 



A free surface boundary layer also exists. Its importance depends on the free 
surface conditions. Consider first the classical case of a clean surface. The 
stress is mainly controlled by the potential velocity field which is of the same 
order as in the main body of the fluid. Because of the small volume 0{6i'^) 
the rate of dissipation in the free surface boundary layer is only 

From the physical point of view it is weaker, since only the zero shear stress 
condition on the free surface is required. 

Another extreme case is when the free surface is heavily contaminated, for 
example, by oil slicks. The stress in the free surface boundary layer can then 
be as great as in the boundary layer near a solid wall. In the present study we 
do not treat such extreme situations and the surface contamination is assumed 
to be absent. 

The previous scalings suggest the following diagram which represents the hi- 
erarchy of dissipative terms: 




It is clear that the largest energy dissipation takes place inside the bottom 
boundary layer. We take into account only two first phenomena from this 
diagram. Consequently, all dissipative terms of order 0{^'^) and higher will 
be neglected. 

Remark 1 In laboratory experiments, surface waves are confined by side walls 
as well. The rate of attenuation due to the side-wall boundary layers was com- 
puted in [ML73]. For simplicity, in the present study we consider an unbounded 
in horizontal coordinates domain (see Figure 1). 
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3 Derivation 



Consider the linearized 3D incompressible NS equations describing free-surface 
flows in a fluid layer of uniform depth h: 

= —Vp + uAv + l V-t7=0, (3) 
ot p 

with V the velocity vector, p the pressure, p the fluid density and g the acceler- 
ation due to gravity. We represent v = {u, v, w) in the form of the Helmholtz- 
Leray decomposition: 

W = V0 + VX^/^, ^/^= (V'l,^/'2,^3). (4) 

After substitution of the decomposition (4) into (3), one notices that the 
equations are verified provided that the functions cf) and satisfy the following 
equations: 

A0 = 0, 0i H ^ gz = 0, — = uAijj. 

p ot 



Next we discuss the boundary conditions. We assume that the velocity field 
satisfies the conventional no-slip condition at the bottom v\_,^_i^ = 0, while at 
the free surface we have the usual kinematic condition, which can be stated 
as 

rit + V ■ V?7 = w. 
After linearization it becomes simply rji = w. 

Dynamic condition states that the forces must be equal on both sides of the 
free surface: 

[£ ■ n] = —(p — po)n + z - n = at 2; = r]{x, t) , 

where a_ is the stress tensor, [/] denotes the jump of a function / across the 
free surface, n is the normal to the free surface and z the viscous part of the 
stress tensor g_. 

The basic idea consists in expressing the vortical part of the velocity field 
V X ■?/; in terms of the velocity potential and the free surface elevation 77 
using differential or pseudodifferential operators. In this section we just show 
final results while the details of computation can be found in [DD07b, Dut07]. 

Let us begin by the free-surface kinematic condition 



Using the absence of tangential stresses on the free surface, one can replace 
the rotational part in the kinematic boundary condition: 



r]t = 4>z + 2uAr]. 



In order to account for the presence of viscous stresses, we have to modify 
the dynamic free-surface condition as well. This is done using the balance of 
normal stresses at the free surface: 

a^z = Uat z = 0^p — Pq = zpu—- = 2pu 



dz ydz"^ dxdz dydz 

In [DD07b] it is shown that £g - = 0{u^), so Bernoulh's equation 
becomes 

0^ + ^^ + 2z/0,, + O(z/i) = 0. (5) 

Since we only consider weak dissipation [u ~ 10^^ — 10~^ m^/s), we neglect 
terms of order o(z/). 

The second step in our derivation consists in introducing a boundary layer 
correction at the bottom. In the present study we assume the boundary layer 
to be laminar. The question of turbulent boundary layer was investigated in 
[Liu06] . Hence, the bottom boundary condition becomes 



90 

dz 



^ f V|0oL=_h , } VOzz 



TT J \Jt — T \ TT J \Jt — r 



dr = -J- '""r'— ^ dr. (6) 







One recognizes on the right-hand side a half-order integral operator. The last 
equation shows that the effect of the diffusion process in the boundary layer 
is not instantaneous. The result is cumulative and it is weighted by (t — r)~2 
in favour of the present time. 

Summarizing the developments made above and generalizing our equations 
by including nonlinear terms, we obtain a set of viscous potential free-surface 
flow equations: 

A0 = O, {x,z) eQ = R'^ X [-h,r]] (7) 

?7t V?7 ■ V</) = 0^ + 2z/A?7, z = ri (8) 

0t + i |V0|' + ^r/=-2z/0,„ z = r] (9) 

, dr, z = -h. (10) 

At the present stage, the addition of nonlinear terms is rather a conjecture. 
However, a recent study by Liu et al. [LPC07] suggests that this conjecture 



t 
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is rather true. The authors investigated the importance of nonhnearity in 
the bottom boundary layer for a sohtary wave solution. They came to the 
conclusion that "the nonlinear effects are not very significant". 

Using this weakly damped potential flow formulation and described in previous 
works procedure [DD07a, Section 4] of Boussinesq equations derivation, one 
can derive the following system of equations with horizontal velocity ue defined 
at the depth zq = -Oh, < 6 < 1: 

rj, + V-{{h + r^)ug) + h'(^^-9 + ^^V\V- ue) = 2uAr^ + dr, 

(11) 

Uet + ^V|Me|' + gVr] - h^O - V(V ■ uet) = 2uAue. (12) 

For simplicity, in this study we present governing equations only on the flat 
bottom, but generalization can be done for general bathymetry. Khabakhpa- 
shev [Kha87], Liu & Orfila [LO04] also derived a similar set of Boussinesq 
equations in terms of depth-averaged velocity. Our equations (11) - (12) have 
local dissipative terms in addition and they are formulated in terms of the ve- 
locity variable defined at arbitrary water level that is beneficial for dispersion 
relation properties. 



3. 1 Total energy decay in visco-potential flow 



Let us consider a fluid layer infinite in horizontal coordinates x = {x,y), 
bounded below by the flat bottom z = —h and above by the free surface. 
Total energy of water waves is given by the following formula: 

S = jj j -\\/<P\^dzdx + ^-jj rfdx. 

E2 -h IR2 

We would like to mention here that Zakharov showed [Zak68] this expression 
to be the Hamiltonian for classical water wave problem with suitable choice 
of canonical variables: rj and ip := 0(5;, z = ri,t). 

In this section we are interested in the evolution of the total energy S with 
time. This question is investigated by computing the derivative ^ with respect 
to time t. Obviously, when one considers the classical potential free surface 
flow formulation [Lam32], we have ^ = 0, since no dissipation is introduced 
into the governing equations. We performed the computation of ^ in the 
framework of the visco-potential formulation and give here the final result: 
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R2 M2 

In this identity, the first term on the right hand side comes from the boundary 
layer and is predominant in the energy decay since its magnitude scales with 
0{y/i'). The second term has its origins in free surface boundary conditions. 
Its magnitude is C(z/), thus it has less important impact on the energy balance. 
This topic will be investigated further in future studies. 

3.2 Dissipative KdV equation 

In this section we derive a viscous Korteweg-de Vries (KdV) equation from just 
obtained Boussinesq equations (11), (12). Since KdV-type equations model 
only unidirectional wave propagation, our attention is naturally restricted to 
ID case. In order to perform asymptotic computations, all the equations have 
to be switched to nondimensional variables as it is explained in [DD07a, Sec- 
tion 2]. We find the velocity variable u in this form: 

U0 = 7] + eP + fi^Q + . . . 

where e, /i are nonlinearity and dispersion parameters respectively (see [DD07a] 
for their definition), P and Q are unknown at the present moment. Using the 
methods similar to those used in [DD07a, Section 6.1], one can easily show 
that 




This result immediately yields the following asymptotic representation of the 
velocity field 

ug = r] - -er]^ + fi^(e - -- —^^1]^^ + ... (13) 

Substituting the last formula (13) into equation (11) and switching again to 
dimensional variables, one obtains this viscous KdV-type equation: 

3 1 * 

^ 

This equation will be used in Section 5 to study the damping of linear progres- 
sive waves. Integral damping terms are reasonably well known in the context of 
KdV type equations. Various nonlocal corrections can be found, for example, 
in the following references [BS71, Che68, Mil76, GPT03]. 

A similar nonlocal KdV equation was already derived in [KM75]. They used 
a different scaling in boundary layer which resulted in dissipative term nonlo- 
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cal in space. Later, Matsuuchi [Mat 76] performed a comparison of numerical 
computations with their model equation against laboratory data. They showed 
that their model does not reproduce well the phase shift: 

... it may be concluded that our modified K-dV equation can describe the 
observed wave behaviours except the fact that the phase shift obtained by 
the calculations is not confirmed by their experiments. 

Excellent performance of the model (14) with respect to experiments was 
shown in [LSVO06]. 



4 Dispersion relation of complete and Boussinesq nonlocal equa- 
tions 

Interesting information about the governing equations can be obtained from 
the linear dispersion relation analysis. In this section we are going to analyse 
the new set of equations (7)-(10) for the complete water wave problem and 
the corresponding long wave asymptotic limit (11), (12). 

To simplify the computations, we consider the two-dimensional problem. The 
generalization to higher dimensions is straightforward and is performed by 
replacing the wavenumber k by its modulus in vectorial case. Traditionally 
the governing equations are linearized and the bottom is assumed to be fiat. 
The last hypothesis is made throughout this study. After all these simplifica- 
tions the new set of equations becomes 

(j)xx + (pzz = 0, (x,z) gM X [-/i,0], (15) 
r]t = (j)z + 2ur]^^, z = 0, (16) 
(l)t + gv + 2u(f),, = 0, z = 0, (17) 

^ 

The next classical step consists in finding solutions of the special form 

0(x, z, t) = if{z)e'^''^-'''^ , r]{x, t) = r/oe*^'^^-^*) . (19) 

From continuity equation (15) we can determine the structure of the function 
(p{z): 

Altogether we have three unknown constants^ C = (Ci,C2,?7o) and three 
boundary conditions (16)-(18) which can be viewed as a linear system with 

^ Since the present problem is linear, we have effectively only two degrees of freedom 
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respect to C: 

MC = 0. 

The matrix M has the following elements 



M 



k 

iuj — 2uk'^ 



iuj — 2uk'^ 



iuj — 2uk'^ 



V 



kF{t,uj] 



^kh 



kF{t,uj] 



where the function F{t, u) is defined in the following way: 



F{t,u;) :-- 



^iul{T){t-T) 



dr. 



(20) 



(21) 



In order to have nontrivial solutions of (15)-(18), the determinant of the sys- 
tem (20) has to be equal to zero det M = 0. It gives us a relation between u 
and wavenumber k. This relation is called the linear dispersion relation: 



D{uj, k) := {iu - 2uk ) + gkiaxi\i{kh) 



kF{t, uj) ({iuj - 2vk'^f tanh(A;/i) + gk^ = 0. (22) 



A similar procedure can be followed for Boussinesq equations (11), (12). We 
do not give here the details of the computations but only the final result: 

Dh{uj, k) := {iuj - 2uk'^f + h{kh)Huj{iuj - 2uk'^) 

+ ghk'^{l - a{khf) - gk^F{t, to) = 0, (23) 

where we introduced the following notation: a := ^ " ^ + |) ^ •= ^(l ^ f ) • 

Unfortunately, the relations D{u,k) = and Dh{uj,k) = cannot be solved 
analytically to give an explicit dependence oi u on k as for the classical wa- 
ter wave problem. Consequently, we apply a quadrature formula to discretize 
the nonlocal term F(t„,a;), where tn denotes the discrete time variable. The 
resulting algebraic equation with respect to uj(tn] k) is solved analytically. 

Remark 2 Contrary to the classical water wave problem and, by consequence, 
standard Boussinesq equations (their dispersion relation can be found in [DD07a, 
Section 3.2], for example) where the dispersion relation does not depend on 
time 

uj'^ - gkiim\i{kh) = (24) 

here we have additionally the dependence of uj{k; t) on time t as a parame- 
ter. It is a consequence of the presence of the nonlocal term in time in the 

but it is not important for our purposes. 
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bottom boundary condition (10). Physically it means that the boundary layer 
"remembers" the flow history. 

Remark 3 There is one subtle point in the derivation presented above. In 
fact, all computations were performed as if the frequency uj were independent 
from time. Our final result shows that time t appears explicitly in the dis- 
persion relations (22), (23). Developments made above make sense under the 
assumption of slow variation of uj with time t. This statement can be writ- 
ten in mathematical form ^ <^ 1. It is rather a conjecture here and will 
be examined in future studies. We had to make this assumption in order to 
avoid complicated integro- differential equations and, consequently, simplify the 
analysis. 

4.1 Analytical limit for infinite time 

In the previous section we showed that the dispersion relation of our visco- 
potential formulation is time dependent. It is natural to ask what happens 
when time evolves. Here we compute the limiting state of the dispersion curves 
(22), (23) as t ^ +00. Namely, we will take this limit in equations (22), (23) 
assuming, of course, that it exists 



Time t comes in dispersion relations through the argument of the function 
F{t,u) defined in (21). Its limit can be easily computed to give 



3 uJoc{k) := lim uJt{k). 



t— >+oo 



lim F{t,uj) 




00 



+ 00 

/5* 




Now we are ready to write down the final results: 



D{u;oo, k) := {iu, 



'00 



2vk'^f + gktwh{kh) 




00 



2vk'^f ta.nh{kh) + gk^ = 0, 



Dh{uJoo, k) := {iuj, 



00 



2h'k'^y + b{khyiuJoc{i^^i 



00 



2iyk'^ 
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In order to solve numerically nonlinear equation D{uJooi k) = (or Dh{uJoo, k) = 
when one is interested in Boussinesq equations) with respect to uJoo, we apply 
a Newton-type method. The iterations converge very quickly since we use 
analytical expressions for the Jacobian (^^) correspondingly). Derivative 
computation is straightforward. Limiting dispersion curves are plotted (see 
Figure 7) and discussed in the next section. 



4-2 Discussion 



Numerical snapshots of the nonclassical dispersion relation ^ at different times 
for complete and Boussinesq equations are given on Figures (2)-(6). The value 
of the eddy viscosity u is taken from Table 1 and we consider a one meter depth 
fluid layer {h = 1 m). We will try to make several comments on the obtained 
results. 



Remark 4 Recall that these snapshots were obtained under the assumption 
thatuo is slowly varying in time. The validity of this approximation is examined 
and discussed in [Dut08]. 

Just at the beginning (when t = 0), there is no effect of the nonlocal term. 
This is why on Figure 2 new and classical ^ curves are superimposed. With no 
surprise, the phase velocity of Boussinesq equations represents well only long 
waves limit (let us say up to kh ^ 2). When time evolves, we can see that the 
main effect of nonlocal term consists in slowing down long waves (see Figures 
3-5). Namely, in the vicinity of kh = the real part of the phase velocity 
is slightly smaller with respect to the classical formulation. From physical 
point of view this situation is comprehensible since only long waves "feel" the 
bottom and, by consequence, are affected by bottom boundary layer. On the 
other hand, the imaginary part of the phase velocity is responsible for the wave 
amplitude attenuation. The minimum of lmCp(/c) in the region of long waves 
indicates that there is a "preferred" wavelength which is attenuated the most. 
In the range of short waves the imaginary part is monotonically decreasing. In 
practice it means that high-frequency components are damped by the model. 
This property can be advantageous in numerics, for example. On Figure 6 we 
depicted the real part of Cp{k) with zoom made on long and moderate waves. 
The reader can see that nonlocal full and Boussinesq equations have similar 
behaviour in the vicinity of kh = 0. 



To be precise, we plot real and imaginary parts of the dimensionless phase velocity 
which is defined as Cr,(k) := — ^^4r^- 

^ In this section expression "classical" refers to complete water wave problem or its 
dispersion relation correspondingly. 
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Figure 2. Real and imaginary part of dispersion curve at t = 0. At the beginning, the 
nonlocal term has no effect. Thus, the real parts of the classical water wave problem 
and new set of equations are exactly superimposed on this figure. The imaginary 
part represents only local dissipation at this stage. 

Now let us discuss the limiting state of phase velocity curves as t ^ +00. It is 
depicted on Figure 7. One can see singular behaviour in the vicinity of zero. 
This situation is completely normal since very long waves are highly affected 
by bottom boundary layer. 

We would like to comment more on the behaviour of curves on Figure 7 since it 
is not easy to distinguish them with graphical resolution. We will concentrate 
on the upper image because everything is clear below with imaginary part. In 
the vicinity oi kh = we have the superposition of nonlocal models (complete 
set of equations (7) - (10) and Boussinesq equations (11) - (12)). When we 
gradually move to short waves, we have the superposition of complete clas- 
sical and nonlocal (7) - (10) water wave problems. Meanwhile, Boussinesq 
system gives slightly different phase velocity for kh > 3. This is comprehen- 
sible since we cannot simplify considerably the problem and have uniformly 
good approximation everywhere. Various Boussinesq systems are designed to 
reproduce the behaviour of long waves. 



5 Attenuation of linear progressive waves 



In this Section we investigate the damping rate of linear progressive waves. 
Thus, the first step consists in linearizing dissipative KdV equation (14) to 
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Figure 3. Phase velocity at t = 1. Boundary layer effects start to be visible: the 
real part of the velocity slightly drops down and the straight lines of the imaginary 
part are deformed by the nonlocal term. Within graphical accuracy, the classical 
and their nonlocal counterparts are superimposed. 




Figure 4. Phase velocity at i = 2. Nonlocal term slows down long waves since the 
real part of the phase velocity decreases. 

obtain the following nonlocal Airy equation: 

^ 

In other words, we can say that we restrict our attention only to small ampli- 
tude waves. Now we make the next assumption. We look for a particular form 
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Figure 5. Phase velocity at t = 4. On Figure 6 we plot a zoom on long and moderate 
waves. 
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(a) Zoom on long waves. Classical and 
corresponding Bousslnesq equations are 
almost superimposed in this region. 



0.7 
0.68 
0.66 
0.64 
' 0.62 

0.6 
0.58 



Nonlocal equations 

- Classical equations 

- Nonlocal Bousslnesq 
■ Classical Bousslnesq 



2.4 2.6 
kh 



(b) Zoom on moderate wavelengths. In 
this region classical and nonlocal com- 
plete equations are almost superim- 
posed. In the same time one can notice 
a little difference in Bousslnesq models. 



Figure 6. Real part of the phase velocity at t = A. 



of the solutions: 

7]{x,t)=A{t)e''^, i = x-^t. (26) 

where k is the wavenumber and A{t) is called the complex amplitude, since 
\ri{x,t)\ = \A{t)\. Integro-differential equation governing the temporal evolu- 
tion of A{t) can be easily derived by substituting the special representation 
(26) into linearized KdV equation (25): 
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Figure 7. Limiting steady state of the phase velocity at t ^ -|-oo. 

In our applications we are rather interested in temporal evolution of the ab- 
solute value It is straightforward to derive the governing equation for 
the squared wave amplitude: 

^ + i.e\AitW - ikM / Mt)Ar)-A(t)A(r) ^ ^ 

at \ irh J Wt — T 

^ 

If we denote by Ar{t) and Ai(t) real and imaginary parts of A{t) respectively, 
the last equation can be further simplified: 

^ + 1.kMitW + 2kM I Mt)Mr)-A(t)A(r) ^ ^ 

at y irh J \/t — T 

^ 

Just derived integro-differential equation represents a generalisation to the 
classical equation (1) by Boussinesq [Bou95] and Lamb [Lam32] for the peri- 
odic, linear wave amplitude evolution in a viscous fluid. We recall that novel 
integral term is a direct consequence of the bottom boundary layer modelling. 

Unfortunately, equation (28) cannot be used directly for numerical computa- 
tions since we need to know the following combination of real and imaginary 
parts Ar(t)Ai{T) — Ai(t)Ar{T) for r G [0,t]. It represents a new and non- 
classical aspect of the present theory. 

Equation (28) allows us to discuss the relative importance of local and nonlo- 
cal dissipative terms for long waves. In fact, when we consider the deep-water 
approximation, only local dissipative terms are present in the governing equa- 
tions [DDZ08]. On the other hand, in shallow waters the integral term is 
predominant. It means that there is an intermediate depth where both dissi- 
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pQjTCLJnctCT 


dcfiTiitiofi 


value 


V 


eddy viscosity 




Q 


gravity acceleration 


9.8 ^ 


h 


water depth 


3600 m 


i 


wavelength 


50 km 


k 


wavenumber 





Table 1 

Values of the parameters used in the numerical computations of the linear progres- 
sive waves amplitude. These values correspond to a typical Indian Ocean tsunami. 

pative terms have equal magnitude. This depth can be estimated when one 
switches to dimensionless form of the equation (28). Comparing the coeffi- 
cients in front of dissipative terms gives the following transcendental equation 
for the "critical" depth h*: 



where uj is the characteristic wave frequency. 

In numerical computations it is advantageous to integrate exactly local terms 
in equation (27). It is done by making the following change of variables: 

^(t) = e-2-^'*e^v^('='^)'*^(t). 
One can easily show that new function A{t) satisfies the following equation: 

— = ik\ — / , A[t) dT 

dt \ nhJ \ft^ ^ ' 

^ 



On Figure 8 we plot a solution of integro-differential equation (27). All pa- 
rameters related to this case are given in Table 1. These values were chosen to 
simulate a typical tsunami in Indian Ocean [DD06]. We have to say that the 
wave amplitude damping is entirely due to the dissipation in boundary layer 
since local terms are unimportant for sufficiently long waves. It means that 
classical formula (2) gives almost constant value ao of the amplitude a(t) on 
the time scale of several hours, since the factor Ivk'^ is of order ^ 10"^^ for 
parameters given in Table 1. 



6 Solitary wave propagation 



In this section we would like to show the effect of nonlocal term on the solitary 
wave attenuation. For simplicity, we will consider wave propagation in a ID 
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1.02 



Amplitude of linear progressive waves 




Figure 8. Amplitude of linear progressive waves as a function of time. Values of all 
parameters are given in Table 1. 

channel. The question of the bottom shear stress effect on the solitary wave 
propagation was considered for the first time in [Keu48]. 

For numerical integration of equations (11), (12) we use the same Fourier- 
type spectral method that was described in [DD07a, Section 5]. Obviously 
this method has to be slightly adapted because of the presence of nonlocal 
in time term. We have to say that this term necessitates the storage of V ■ 
M*^"^ at previous time steps. Hence, long time computations can be memory 
consuming. 

The values of all dimensionless parameters are given in Table 2. Dimensionless 
viscosity v is related to other physical parameters in the following way: 



V 



V 



gh 



where i> is kinematic viscosity, i is the characteristic wavelength and h is the 
typical depth. 

Remark 5 From numerical point of view, the integral term J '^^^(^'~^'^) dj- 

^ 

can pose some problems in the vicinity of the upper limit r = t. Probably the 
best way to deal with this problem is to separate the integral in two parts: 



I 



dr 



dr 



,(f, -h,T) 



dr, 6 > 0. 



The first integral can be computed in a usual way without any special care. 
Then one applies to the second integral a special class of Gauss quadrature 
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pCtTCLUlCtCT 


dcfifiition 


valuG 


£ 


nonlinGarity 


u.uz 




dispersion 


U.UD 


u 


eddy viscosity 


0.001 


C 


soliton velocity 


1.02 


9 


zq = —Oh 


1 - \/5/5 




soliton center at t = 


-1.5 



Table 2 

Values of the dimensionless parameters used in the numerical computations. 
formulas with weighting function (t — r)^^ . 

But there is another well-known trick that we describe here. This technique 
can be implemented in simpler way. We rewrite our integral in the following 
way 

j (pzzjx, -h,T) f <Pzzix, -h,t) f <Pzzix, -h, t) - (pzzjx, -h, t) 

J ^/t^ ^ ~ J ^/t^ ^ J Vt^ ^' 

^ 

The first integral in the right hand side can be evaluated analytically while the 
second one does not contain any singularity under the assumption of differen- 
tiability of r ^ (pzzix, —h,r) at r = t: 

J M^^H dr = 2v^0..(f , -h, t) + J <I^US^-h,r)_^M-h,t) 
^ ^ 



Remark 6 What is the value ofu to be taken in numerical simulations? There 
is surprisingly little published information of this subject. What is clear is that 
the molecular diffusion is too small to model true viscous damping and one 
should rather consider the eddy viscosity parameter. Some interesting infor- 
mation on this subject can be found in [TSL07]: 

We have spent a considerable amount of time and effort seeking further 
published information on viscous effects on ship waves. 

This sentence confirms our apprehension. The authors of this work came to 
the following conclusion 

...we reiterate that a viscosity of about v = 0.005 m^/s gave reason- 
able agreement with longitudinal cut results (including apparent damping 
of transverse waves). 
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In another famous paper [BPS81] one can find: 

. . . Such a decay rate leads to a value for n in (M* ) of 0.014- ■ ■ 

In their work fi is the coefficient in front of the dissipative term in BBM 
equation: 

3 1 

Vt + 'nx + T^vvx - i^Vxx - -^Vxxt = 0. 

It is important to underline that this equation is written in dimensionless 
variables. Thus, the value reported in that study has to be rescaled with respect 
to other physical parameters. 

Liu and Orfilla [LO04] report the value of eddy viscosity v = 0.001 m? /s. 

If we summarize all the remarks made above, we can conclude that the value 
of the order 10""^ - 10~^ m? /s should give reasonable results. 

6.1 Approximate solitary wave solution 

In order to provide an initial condition for equations (11), (12), we are going 
to obtain an approximate solitary wave solution for nondissipative ID version 
of these equations over the flat bottom: 

r]t+ ({1 + eri)u)^ + j/i— - + - \u^xx = ^, 
ut + ilx + ^iu^)x - /i^6'(^l - ^^u^xt = 0- 

Then, we apply the same approach as in Section 3.2 or in [DD07a, Section 
6.1]. We do not provide the computations here since they are simple and can 
be done without any difficulties. The final result is the following: 

r]{x, t) = — ^sech ^ ^„ —{x + Xq — at)) 

and the velocity is given by formula (13). In the numerical results presented 
here, we use ri{x, 0) and u{x, 0) as initial conditions. 

6.2 Discussion 

On Figures 9-11 we present three curves. They depict the free surface eleva- 
tion according to three different formulations. The first corresponds to classical 
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t = 0.500 




X 



Figure 9. Comparison among two dissipative and nondissipative Boussinesq equa- 
tions. Snapshots of the free surface at t = 0.5 



t = 1 .000 




X 



Figure 10. Comparison among two dissipative and nondissipative Boussinesq equa- 
tions. Snapshots of the free surface at t = 1.0 

Boussinesq equations without dissipation. The second one to dissipative sys- 
tem with differential or local terms (for example, vAu in momentum conser- 
vation equation) and the third curve corresponds to the new set of equations 
(11), (12). On Figure 12 we made a zoom on the soliton crest. 

It can be seen that Boussinesq equations with nonlocal term provide stronger 
attenuation of the amplitude. In the same time, as it was shown in the previous 
section, this nonlocal term slightly slows down the solitary wave. 

In order to show explicitly the rate of amplitude attenuation, we plot on Fig- 
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t = 2.000 




X 



Figure 11. Comparison among two dissipative and nondissipative Boussinesq equa- 
tions. Snapshots of the free surface at t = 2.0 

t = 2.000 




X 



Figure 12. Comparison among two dissipative and nondissipative Boussinesq equa- 
tions. Zoom on the sohton crest at t = 2.0 

ure 13 the graph of the following apphcation t sup_.^<^<^ \ri{x,t)\. One can 
see on this plot little oscillations which are of numerical nature. Our numerical 
experiments show that their amplitude decreases when — > oo. This result 
shows one more time that nonlocal model provides stronger damping proper- 
ties. One can have the impression that the amplitude decays linearly but it is 
only an impression because of (2). This behaviour for small t can be explained 
by simple Taylor expansion which is valid when uk'^t -C 1: 

a{t) = aoe-^"'''' = ao{l - 2pkH) + 0(uW). 
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Maximum wave amplitude as function of time 




In the present article we used a novel visco-potential formulation proposed re- 
cently in [LO04, DD07b, Dut07]. This formulation contains the nonlocal term 
in the kinematic bottom boundary condition. This term should be considered 
as a boundary layer correction at the bottom. In modelling viscous effects this 
term plays the main role, since its magnitude is 0{^/l'). Recently, important 
progress was made in numerical computation of the nonlocal term [TL07]. 
Thus, we can hope to have its implementation in operational numerical codes. 
It is interesting to note that the boundary layer effect is not instantaneous 
but rather cumulative. The flow history is weighted by (t — r)"2 in favour of 
the current time t. As pointed out in [LO04], this nonlocal term is essential to 
have an accurate estimation of the bottom shear stress based on the calculated 
wave field above the bed. Then, this information can be used for calculating 
sediment-bedload transport fluxes and, in turn, morphological changes. 



A long wave asymptotic limits (Boussinesq and KdV equations) were derived 
from this new formulation. The dispersion relation was described. Due to the 
presence of special functions, we cannot obtain a simple analytical dependence 
of the frequency u on the wavenumber k as in classical equations. Consequently 
the dispersion curve was obtained numerically by a Newton-type method. We 
made a comparison between the phase velocity of the complete visco-potential 
problem and the corresponding Boussinesq equations. The dispersion relation 
of new formulation is shown to be time dependent and this property is not 
classical. It comes from the memory effect of the boundary layer. We computed 
analytically the limiting dispersion curve ujoo{k) as t ^ +00. 
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The effect of the nonlocal term on the solitary wave attenuation was inves- 
tigated numerically with a Fourier-type spectral method. We showed that it 
provides much stronger damping (see Figure 13). An equation describing the 
amplitude evolution of linear progressive waves was derived (27). This result 
includes bottom boundary layer correction and generalizes classical formula 
(1) by Boussinesq [Bou95] and Lamb [Lam32]. 

The present study opens new exciting possibilities for future research. Namely, 
we did not consider at all the questions of theoretical justification of visco- 
potential formulation and corresponding long wave models in the spirit of 
works by J. Bona, J.-C. Saut, D. Lannes, T. Colin [BCL05, BCS02]. Other 
directions for future work consist in implementing new visco-potential frame- 
work in operational ocean and nearshore models. 
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